Anisotropic frustrated Heisenberg model on the honeycomb lattice 
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We investigate the ground-state phase diagram of an anisotropic Heisenberg model on the hon- 
eycomb lattice with competing interactions. We use quantum Monte Carlo simulations, as well 
as linear spin-wave and Ising series expansions, to determine the phase boundaries of the ordered 
magnetic phases. We find a region without any classical order in the vicinity of a highly frustrated 
point. Higher-order correlation functions in this region give no signal for long-range valence-bond 
order. The low-energy spectrum is derived via exact diagonalization to check for topological order 
on small-size periodic lattices. 

PACS numbers: 75.10.Jm; 75.30.Kz; 75.40.Mg 
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I. INTRODUCTION 

The honeycomb lattice has received much attention in 
recent years because of its relevance to graphene, whose 
electronic structure gives rise to many unusual features. 1 
However, this two-dimensional bipartite lattice with its 
two-site unit cell was investigated long before it was re- 
alized in a real material. It is particularly interesting for 
quantum-mechanical models of strongly correlated elec- 
trons since its coordination number n = 3 is the lowest 
allowed in a two-dimensional system. 2,3 Hence the influ- 
ence of quantum fluctuations on the ground-state prop- 
erties is expected to be more important than, e.g., in 
the also bipartite square lattice. Recently, it was found 
that the half-filled Hubbard model on the honeycomb 
lattice exhibits a spin liquid state at the border of the 
metal-insulator transition for an intermediate value of 
the on-site repulsion (7. 4,5 Since then the investigation of 
spin models that can be derived perturbatively from the 
Hubbard model, 6 or stated directly inside the insulating 
phase, has yielded many interesting features including 
disordered and valence bond solid phases in the ground- 
state phase diagram. 7-12 Furthermore, in the context of 
iridium compounds, systems have been investigated that 
may be described by frustrated spin models on the hon- 
eycomb lattice and exhibit antiferromagnetic order. 13 

Quantum Monte Carlo (QMC) simulations are rather 
difficult to apply for these spin models due to the sign 
problem that appears in relation to frustrating quan- 
tum spin fluctuations and thermalization problems ac- 
companying the competition between different ground 
states. The latter can be overcome by additional ex- 
change Monte Carlo updates, as was shown in previ- 
ous work for frustrated spin models 1 and in particular 
for an anisotropic J1-J2 Heisenberg model on the square 
lattice. 15 To avoid the sign problem completely and allow 
for the investigation of ground-state properties of reason- 
ably large systems it is necessary to lift the frustration 
for some interactions in such a way that the isotropy of 
the model is lost. 



In the present work we investigate a spin model on 
the honeycomb lattice, including nearest-, next-nearest-, 
and third-nearest-neighbor anisotropic Heisenberg inter- 
actions, a geometry analyzed in previous works for the 
isotropic Heisenberg model. 2,8,10 We start from the limit 
of small quantum fluctuations, where classical antifer- 
romagnetic S z interactions along a preferred direction 
are all antiferromagnetic, leading to frustration. How- 
ever, we choose the quantum fluctuations in the transver- 
sal plane to be ferromagnetic and hence nonfrustrating. 
Such an anisotropic version of the Heisenberg model can 
be interpreted as a system of hard-core bosons 16 with 
nonfrustrating kinetic energy and repulsive interactions. 
In consequence, an application might be realized in opti- 
cal lattices. 17,18 

Applying different techniques, we find the ground-state 
phase diagram of the anisotropic spin model. We analyze 
the boundaries of the ordered phases and explore the re- 
gion between them. Using QMC simulations and exact 
diagonalization (ED), we identify a finite parameter re- 
gion where a disordered ground state cannot be excluded. 
Very similar behavior was observed for the anisotropic J\- 
J2 Heisenberg model on the square lattice. 15 The paper 
is structured as follows. In Sec. II we introduce the spin 
model and its equivalent hard-core boson realization. We 
briefly present in Sec. Ill the different methods that will 
be applied to derive the phase diagram. The results for 
different limiting cases and from the various approaches 
are compared and discussed in the main part of the paper 
in Sec. IV. The possibility of a spin liquid phase, as well 
as advantages and drawbacks of the different methods, is 
discussed in the concluding Sec. V. 



II. MODEL 

We consider a spin model on a periodic honeycomb 
lattice with N = 2 x (L x L) sites. The homogeneous 
anisotropic Heisenberg interactions between spin opera- 
tors S = (S x , S y , S z ) at sites i,j, separated by a distance 
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(b) Schematic phase diagram 

FIG. 1. (a) Classical ground states: Neel and collinear 
order on the honeycomb lattice, (b) Schematic phase di- 
agram for the frustrated anisotropic Heisenberg model on 
the honeycomb lattice. The axes represent the frustration 
V = V2/V1 = V3/V1 and the relative amplitude of the quan- 
tum fluctuations t = t r /V r . Red circles refer to actual transi- 
tion points calculated with QMC simulations (see below). 



called r, are given by 

J r (Sj, SjJ r = J r Si Sj + J r y ^Sj Sj + Sf Sjj i (1) 

where z refers to a preferred spin direction, x and y 
are the transversal directions, and J r = (Jr,Jr ,y ) is the 
anisotropic exchange coupling strength. The Hamilto- 
nian is given by summing such interactions over all spin 
pairs on nearest-neighbor (r = l), next-nearest- neighbor 
(NNN) (r = 2) , and third-nearest-neighbor (r = 3) bonds 
of the honeycomb lattice: 

ff = E E Jrfr,Sj) r - (2) 

r = 1 {i,j)r 

Note that J\ and J3 couple sites of the different sublat- 
tices (say, A and B) of the bipartite honeycomb lattice, 



while the NNN sites coupled via J2 lie in the same sublat- 
tice. We choose the S z interactions to be antiferromag- 
netic, J z > 0, such that in the Ising limit (Jr ,y = 0) Jf 
and Jf favor a Neel configuration, with opposite spin ori- 
entation in each sublattice [see Fig. 1(a), left]. However, 
with Jf > a competing interaction is introduced, gener- 
ating frustration in the model. The sign of the exchange 
terms J*' y and J3 ' v does not matter since a sublattice 
rotation can change both signs simultaneously. However, 
for the NNN exchange a negative (ferromagnetic) sign 
has to be chosen to avoid the sign problem in the QMC 
simulations. For convenience we set all Jr' v < so that 
the fluctuations are ferromagnetic and hence nonfrustrat- 
ing. 

The present anisotropic Heisenberg model, for spin 
S = 1/2, can also be understood as a model of hard- 
core bosons. For convenience we first fix an energy scale 
J, defining dimensionless exchange couplings J*/J = V r 
and Jr' v I J = 2t r . We rewrite the spin model in terms of 
ladder operators = Sf ± iSf , 

3 

^=EE (VrSfSl + tr(S+SJ + SiSfj) (3) 

r=1 {i,j)r 

and map spin operators onto hard-core bosons by (cf., 
e.g., Ref. 16) 

Sf^b\, S-^b^ Sf-^nj-1/2, (4) 

where 

r H = b\b i and (&| f) ) 2 =0, (5) 
so that the Hamiltonian (3) maps to 

3 

H hoson /J = E E { tr ( b l b j + b ] b i) + VrUiTlj) 
r=1 (i,j)r 

- |(Vi + 2V 2 + V 3 )J2^ + f(Vi + 2^2 + V 3 )N . 

i 

(6) 

Here the (negative) t r describe nonfrustrating hopping of 
the hard-core bosons up to third-nearest neighbors and 
the (positive) V r describe repulsion up to the same range. 
In particular we are interested in the zero magnetization 
case, mapping to a half-filled lattice, i.e., (rn) = 1/2. In 
this case the last line in Eq. (6) only yields constant terms 
that will be dropped in all following considerations. 

For the remainder of this work, all three interactions 
will have the same anisotropy ratio Jr' v /Jr = 2t < 0, 
the relative strength of the interactions is set to V2/V1 = 
V3/V1 = V > 0, and the scale is set to J\ = J (i.e., V\ = 1, 
implying t\ — t and t<i — t 3 = Vt). Thus the parameter 
space is reduced to a two-dimensional area, which we will 
investigate for positive V and negative t. 

Ground-state configurations are readily obtained in 
some limiting cases and a schematic phase diagram is an- 
ticipated in Fig. 1(b). In the classical Ising limit (t — > 0) 
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of the spin model in Eq. (3), two different ground states 
occur: for V < 1/2 a Neel state [Fig. 1(a), left], with 
energy 

Emd/J= §(V-l)S 2 iV, (7) 

and for V > 1/2 a collinear state [Fig. 1(a), right], with 
energy 

E coXl/J = \{l-W)S 2 N. (8) 

At the critical point V = 1/2 the two ground states 
compete and the transition temperature is suppressed 
to zero. The degenerate ground-state manifold at this 
point is expected to give rise to interesting phenomena 
for nonzero quantum fluctuations (t < 0) in the direct 
vicinity of the critical point. 

The limit of large values of \t\ can be easily under- 
stood too as the noncompeting ferromagnetic interac- 
tions yield a ferromagnetic correlation in the x-y plane 
[see Fig. 1(b)]: The ground state will spontaneously 
break the rotation symmetry in the x-y plane, still with 
vanishing magnetization in the z axis. The anisotropic 
Heisenberg interactions (1) can also be understood in 
terms of Ising interactions between the, say, S x spin com- 
ponents, regarding the interactions between the other 
spin components S — S y ± iS z as (not negligible) spin 
fluctuations: 

r = 1 {i,j)r 

+ £ £ ^(S+S++H.c.) 

+ EE W^+H.c). (9) 

r=1 (i,3)r 

Then, for large \t\ a ferromagnetic product wave function 
in the S x basis serves as variational ansatz (VA), with 
energy 

E ie „o/J = St(l + 3V)S 2 N . (10) 

In the bosonic language, the limits described above 
match density waves at small \t\ and a superfluid phase 
at large \t\, corresponding to the condensation of ferro- 
magnetic magnons. 16,20,21 Note that Neel and collinear 
phases are only exact ground states for t — or large 
S, i.e., in the classical limit, while the bosonic treat- 
ment is valid for 5 = 1/2 and will be studied at finite t. 
One should then expect that fluctuations will reduce the 
magnetic order. To determine the stability boundaries of 
these phases in the V-t plane [sketched in Fig. 1(b)] and 
to analyze the intermediate region of the phase diagram, 
we apply various methods, which are briefly introduced 
in the following section. 



III. METHODS 

In this section we will briefly summarize the meth- 
ods that we have employed; results will be presented in 
Sec. IV. The reader who is not interested in technical 
details may skip this section. 

The series expansion (SE) method will be applied in 
the limit of small fluctuations to calculate energies and 
estimate the phase boundaries of the antiferromagnetic 
phases. The same holds for the derivation of linear spin 
waves (LSWs), which will also be applied for the vari- 
ational ferromagnetic state. The QMC simulations are 
employed for the whole phase diagram to calculate vari- 
ous order parameters and results for the energies will be 
compared to the other methods. In addition, for a par- 
ticular parameter set, ED will be performed to calculate 
the low-energy spectrum. 



A. Series Expansion 

We have analyzed perturbatively the Hamiltonian (3) 
starting from the Ising limit (t r = 0, r = 1, . . . , 3), around 
classical Neel and collinear phases. Using standard 
(Rayleigh-Schrodinger) perturbation theory 22 on finite 
lattices, we have obtained analytic expressions for the 
ground-state energy up to 0(tf) around the two men- 
tioned classical phases. We have performed these cal- 
culations on finite lattices, large enough to avoid finite- 
size effects at this order, using computer algebra software 
for the implementation. 23 The method is straightforward 
but usually limited to low orders due to the memory con- 
straints imposed by the need of larger lattices to achieve 
higher orders. Although our expansions are only fourth 
order, they have the advantage of providing the coeffi- 
cients of the expansions in a closed analytical form. This 
allowed us to determine a first-order critical line between 
Neel and collinear phases, as discussed in the following 
sections. 

Very recently, Oitmaa and Singh 11 analyzed the 
isotropic Heisenberg model applying a linked-cluster 
formalism 24 and calculated numerically series up to 
eighth order. For such a purpose an Ising model is per- 
turbed with an auxiliary anisotropy parameter that in the 
end has to tend to one to recover the isotropic Hamilto- 
nian. This parameter is directly related to the couplings 
t r in Eq. (3) ; thus their results are useful in analyzing our 
model. We use here their numerical eighth order series 
for two purposes: First, we use the ground-state energy 
results to estimate a range of validity of our series; second 
and more importantly, we use the series provided for or- 
der parameters to determine critical points in the phase 
diagram of our model. 
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B. Linear Spin Waves 

The phase boundaries of the conventional magnetic 
phases can be estimated for large spin S, using a 1/5 
LSW expansion of the model in Eq. (3). To this end one 
selects a classical spin configuration to set local Cartesian 
frames with z% along the classical direction and transver- 
sal components Xi and yi and represents spins in terms 
of Holstein-Primakoff bosons a,; satisfying = ^ij 25 

The spin component along the local classical direction is 
represented by 

Sf = Sj • £i = S — rii, (11) 

where rij = a\ai, while transversal fluctuations Sf- = Sj ■ 
Xi ± iSi ■ iji are represented by 

S+ = (v/25 - m) oi, (12) 

Sf =a\ (v/25-ni) . (13) 

These equations ensure the correct spin component com- 
mutation relations as well as —5 < (Sf) < S. The LSW 
expansion is then obtained by expanding the Hamiltonian 
up to second order in a^; while this procedure usually 
breaks SU(2) invariance, in the present case that symme- 
try is already explicitly broken by the anisotropy. 

The form of the bosonic action depends on the selected 
classical configuration. In the present case we have ex- 
plored the Neel and collinear configurations as well as a 
ferromagnetic phase in the x-y plane. In these three con- 
figurations the local classical directions at different sites 
are either parallel or antiparallel, simplifying the com- 
putations. For classically parallel interacting sites the 
interaction terms in Eq. (3) are expanded as 

V r S 2 - V r S (rii + rij) + 2St r (a\aj + a\a^ , (14) 

while for antiparallel local frames they read 

-V r S 2 + V r S (m + nj) + 2SU [a\a] + a t a^ . (15) 

The quadratic Hamiltonian contains anomalous terms 
in the fluctuations and must be diagonalized by a Bo- 
goliubov transformation. Within each momentum mode, 
the Bogoliubov transformation is possible if the Hamil- 
tonian quadratic form has only positive eigenvalues. We 
searched such regions in the V-t plane and after checking 
that fluctuations do not destroy the classical order we 
adopted the possibility of diagonalizing the Hamiltonian 
as a criterion for stability of the classical phase. Notice 
that the ferromagnetic phase shows a Goldstone mode 
at zero momentum, related to the explicit selection of a 
classical direction in the x-y plane; in this phase, stability 
is achieved if the remaining modes are positive. 

Within each phase, we computed the ground-state en- 
ergy and the average magnetization along the classical di- 
rections as order parameters. In the regions where more 
than one phase is stable, we select the one with lowest 
ground-state energy. The results are shown in Sec. IV. 



C. Quantum Monte Carlo 

The simulation of frustrated spin models using QMC 
simulations is usually limited by the sign problem that 
occurs for competing quantum fluctuations. However, for 
the present model we choose the interactions along the 
S xy direction of the spin operators to be ferromagnetic 
and hence have negative amplitudes t r for the quantum 
fluctuations. This maps onto bosonic hopping integrals 
[see Eq. (6)], which yield no sign problem in QMC sim- 
ulations. Thus we can simulate the model for all pa- 
rameters. Nevertheless, the remaining frustration in the 
S z direction of the spin interactions yields a competition 
between different ground states for V « 1/2. Exactly 
at the critical point the classical ground state shows a 
large degeneracy of linear order ~ 2 2 L and this results in 
freezing and thermalization problems for simulations of 
the quantum model at low temperatures. 

The QMC simulations were performed at finite temper- 
atures using the implementation of the stochastic series 
expansion 26 of the ALPS project. 27-29 The use of directed 
loops 30 in the update step ensures a reliable scan of the 
whole phase space even for the present model where sev- 
eral different operators can act on the same site of the lat- 
tice. To overcome the freezing problems that occur due to 
the degenerate ground-state manifold at V = 1/2 we used 
an additional exchange Monte-Carlo update ' 31 ' 32 that 
allows for a better thermalization of the simulation. Since 
in particular in the vicinity of the critical point V « 1/2 
low temperatures are necessary to gain insight into the 
ground-state properties of the system, simulations were 
run in parallel on large-scale computer clusters. 

D. Exact Diagonalization 

To gain further insight into the spectrum of the model 
we also applied exact diagonalization techniques. This is 
done for finite lattices up to N = 34 sites and only for cer- 
tain parameters where SE and LSWs are not applicable. 
In particular we performed Lanczos iterations to calcu- 
late the lowest eigenvalues of the system using an imple- 
mentation by Jorg Schulenburg. 33 The finite-size analysis 
of the behavior of the lowest excited states allows for a 
characterization of the underlying ground state. 

IV. RESULTS 

In this section we present and compare results obtained 
by the methods described above for three different cases 
of interest. For small fluctuations the stability of the 
antiferromagnetic states will be analyzed by means of 
SE, LSWs and QMC simulations. In the case of large 
fluctuations the emergence of long-range ferromagnetic 
correlations in the x-y plane is tested by means of LSWs 
and QMC simulations through the spin stiffness of the 
S z component. An intermediate region in the vicinity of 



5 



-0.2 
-0.3 
-0.4 : 
-0.5 





□ QMC 

LSW 

SE (4th order) 

X SE (8th order, Ref. 11) 




FIG. 2. Comparison of energies calculated from the different 
methods introduced in Sec. Ill, at small quantum fluctuations 
t = -0.10 (top) and t = -0.05 (bottom). For the SE two dif- 
ferent calculations up to fourth order (this work) and up to 
eighth order (Ref. 11) are shown. In the direct vicinity of the 
critical point V — 1/2 both expansions become rather unre- 
liable due to an increasing number of divergences. The LSW 
results underestimate the influence of quantum fluctuations, 
as also observed for larger \t\. 



the critical point V = 1/2 will be investigated via QMC 
calculations of higher-order correlation functions and in- 
terpretation of the low-energy spectrum from ED. 
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FIG. 3. Ground-state phase diagram for the anisotropic 
Heisenberg model on the honeycomb lattice by means of QMC 
simulations (red circles and interpolation by red dashed lines). 
The Ising ground states survive separated by a direct first- 
order transition for small fluctuations \t\. Only for values of 
\t\ > 0.175(25) can ferromagnetic order in the x-y plane be de- 
tected. The solid blue line represents the first-order transition 
line between the antiferromagnetic states [Eq. (16)], deter- 
mined from fourth order SE. Magenta crosses represent phase 
boundaries determined by the condition of vanishing order pa- 
rameters, provided as eighth order SE (Ref. 11) (see the text 
for more details). Inset: The LSWs (green dash-dotted lines) 
yield phase transition lines due to stability arguments and a 
comparison of energies between antiferromagnetic states. The 
overall scenario is very similar. 



The Ising limit is given by setting all quantum fluctu- 
ations t r — and exhibits two antiferromagnetic ground 
states for Va = V3 that were described in Sec. II. For small 
values of \t r \ the quantum- mechanical ground states are 
expected to consist of those classical states plus some 
quantum fluctuations that reduce the overall energy and 
order parameters. 

A comparison of the overall energies from SE, LSWs, 
and QMC simulation is given in Fig. 2 for fluctuations 
governed by t = —0.05 and —0.10. The agreement for 
small V < 0.45 and large V > 0.7 is very good. Only in 
the intermediate regime can discrepancies be observed, 
which will be discussed below. 

From our fourth-order SE calculations we evaluated 
ground-state energies for the Neel (V < 1/2) and collinear 
states (V > 1/2). Series for these two states are given in 
the Supplemental Material. 34 These results agree with 
the recent work by Oitmaa and Singh 11 up to the given 
order. 

The LSW approximations yield comparable results for 
the energies that are shown in Fig. 2 for t = —0.05 and 
—0.10, computed on a lattice with 2x 10 4 sites. The intrin- 
sic breakdown of the method for a particular starting con- 
figuration as described in Sec. Ill B provides an estimate 
of the upper phase boundary for both antiferromagnetic 
states: The Neel configuration is stable for —t < j-f^y 



and the collinear one is stable for — t < — fe^rr ■ Both an- 
tiferromagnetic configurations support LSW fluctuations 
for V w 1/2, where the phase boundary is estimated by 
energy comparison. The three approximate LSW bound- 
aries are plotted in the inset of Fig. 3. 

Our most accurate results were obtained by extensive 
QMC simulations for the energies and magnetic order 
parameters. To identify the regions with different mag- 
netic orders we calculate the structure factors for the 
Neel and collinear order. The Cartesian wave vector is 
given by k = (0, 0) and antiparallel spins on the A and 
B sites of the unit cell for the Neel state, i.e., each sub- 
lattice is ferromagnetically ordered but they are aligned 
antiparallel to each other. The collinear state is sixfold 
degenerate with three wave vectors: k = -^=(-\/3, 1) and 

k = ^(V3, -1) with A and B parallel and k = ^(0, 1) 
with A and B antiparallel. Additionally all spins can be 
flipped in the ordered states, which gives a twofold degen- 
eracy. An example of the behavior of the order parameter 
as a function of temperature is shown in Fig. 4(a) for the 
Neel state. 

As expected from SE for small \t\, we obtain a direct 
transition between the two antiferromagnetic states that 
is probably of first order since both states exhibit dif- 
ferent symmetries. From QMC simulations we observe 
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FIG. 4. Evolution of QMC results for energies and magnetic 
order parameters shown for decreasing temperatures in two 
different points of the phase diagram, exhibiting (a) Neel or- 
der (NO) and (b) ferromagnetic order (FO), respectively. 



that this transition line splits into two for a small value 
of 0.025 < \t\ < 0.05 and a new ground state emerges in 
between. 



B. Ferromagnetic limit 

In the opposite limiting case with large \t\ the situa- 
tion is not as simple as for the Ising limit. Linear spin 
waves can be expanded around the product state that 
yields the energy given in Eq. (10) and again the break- 
down will yield an estimation of a lower phase boundary. 
However, the coupling strength of the quantum fluctua- 
tions, scaling with \t\, is not small compared to the Ising 
exchange and the position of this phase boundary is not 
very reliable. In the inset of Fig. 3 we show the ferro- 
magnetic phase boundary obtained by LSWs in a large 
lattice, well fitted by -t > ■ Then LSWs cannot 

be applied in the unidentified region in the inset of Fig. 3. 
The energies obtained by this calculation are compared 
below to the results of QMC simulations. 

The appropriate order parameter in the QMC simu- 



0.2 



^ t = - 0.05 



- -0.2 



-o.i-Er' 0--OFO \ □ QMC 
A: A no v 

□■ -dCO 



a a energy / Yn 

A-ANO • __ SE (4th order) , W 

/ 




-0.22 3 



-0.24 53 



0.4 



(a) fixed t = —0.05 and varying frustration V, L = 12 



V = 0.45 




(b) fixed V = 0.45 and varying fluctuations t, L = 12 

FIG. 5. Energy and order parameters (FO, ferromagnetic in- 
plane order; NO, Neel order; CO, collinear order) shown for 

(a) a horizontal cut through the phase diagram (Fig. 3) at 
t = —0.05 (classical order is absent for 0.48 < V < 0.54) and 

(b) a vertical cut at V = 0.45 (no finite order parameter for 
0.075 < -t < 0.15). 



lations for the ferromagnetic state is given by the spin 
stiffness, which can be estimated on behalf of the winding 
number within the QMC algorithm. 16,35 As an example, 
the convergence of energies and this order parameter is 
shown in Fig. 4(b). Careful calculations of this observ- 
able for the not antiferromagnetically ordered regions of 
the phase diagram show only a nonvanishing signal for 
\t\ > 0.15. Figure 5 shows two different parameter scans. 
In Fig. 5(a) order parameters and energies from QMC 
simulations and SE are shown for t = —0.05 and varying 
frustration V from Neel to collinear behavior: A finite re- 
gion without any magnetic order is identified, which also 
explains the discrepancy of the energies shown here and 
in Fig. 2. In Fig. 5(b) a similar scan is presented, here 
for fixed V = 0.45 and t varying from Neel to ferromag- 
netic behavior, where the QMC energy is also compared 
with the classical variational ansatz [Eq. (10)] and LSW 
calculations. The agreement with the LSW energy is re- 
markably good. 

As a result from the QMC, SE and LSW calculations 
we show a phase diagram for the V-t parameter space in 
Fig. 3. The three methods uncover a finite region with- 
out magnetic order around the critical point V = 1/2 and 
finite fluctuations — t, which will be analyzed in the fol- 
lowing section. The Ising ground states survive separated 
by a direct first-order transition for small fluctuations \t\. 
Only for values of \t\ > 0.175(25) is ferromagnetic order in 



the x-y plane detected by means of QMC simulations. A 
direct comparison of energies for V rs 1/2 yields a phase 
transition line between the two antiferromagnetic states. 
We give here the approximate result as a function V(t), 
obtained by equating our fourth order SE of energies for 
both Neel and collinear states and expanding up to linear 
order in V, around V = 1/2: 



„• Q„ 



V(t) 



1278676 1_ - 69750 r + 6300 1_ - 225 
2665577 1 4 - 148050 1 3 + 13950 1 2 - 450 



(16) 



This function is shown as solid blue line in the ground- 
state phase diagram in Fig. 3. 

As we have mentioned, Oitmaa and Singh 11 have pre- 
sented SE data for the magnetic order parameters, i.e., 
Neel and collinear magnetization, that can be directly 
translated to our model. We used these data to deter- 
mine the upper phase boundaries of the antiferromag- 
netic phases by detecting the point in V-t space where 
order parameters (more precisely, their Pade approxi- 
mants) vanish. The result is also shown in the phase 
diagram (Fig. 3, magenta crosses). The corresponding er- 
ror bars are confidence limits obtained by considering the 
dispersion of predicted critical points for different Pade 
approximants. Note that the antiferromagnetic phase 
boundaries estimated by SE are generally slightly above 
the numerical QMC results. A possible explanation is the 
presence of first-order transitions at these phase bound- 
aries. 

The LSWs (green dash-dotted lines in the inset of 
Fig. 3) yield phase transition lines as described in 
Sec. IIIB. From the Monte Carlo simulations, the phase 
diagram is obtained by identifying the regions where the 
different order parameters show finite signals (red circles 
in Fig. 3). As shown in Fig. 5, this leaves a finite region 
with a ground state that shows no magnetic order. Since 
the energies calculated by LSWs were slightly above the 
values of QMC simulations and SE, for increasing — t it 
is not surprising that the phase boundaries are shifted to 
higher values as well. 



C. Intermediate regime 

For the intermediate regime we find similar behavior 
as for the square lattice, 15 i.e., only small signals ap- 
pear in the spin stiffness for small lattices and interme- 
diate temperatures that scale to zero for larger lattices. 
The estimation of the second critical value for t(V), at 
which ferromagnetic order arises, is rather difficult due 
to the two-dimensional parameter space and the finite- 
size problems in the ferromagnetic order parameter. In 
Fig. 3 we show a rough estimate for the phase transition 
line separating the disordered state and the ferromag- 
netic phase as a red dashed line. 

So far we only checked the disordered region for classi- 
cal magnetic order. However, frustrated quantum models 
are known to exhibit in certain cases more exotic quan- 
tum ordered patterns such as dimer phases with long- 
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(a) Neel state at V = 0.2 and t = -0.1 (T = 0.05 J) 



0.31832 
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0.00 
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(b) ferromagnetic state at V = 0.45 and t - 
,«- as O. „« "'« 



-0.5 (T = 0.1 J) 



0.09111 



,0°°°'0„ 



-0.00690 
0.25 



(c) disordered state at V = 0.45 and t = -0.1 (T = 0.05 J) 

FIG. 6. QMC results for correlation functions on an N — 288 
lattice. Each bond (and site) represents the strength of the 
correlation of a dimer (or spin) at the corresponding distance 
to the top left dimer (or spin) in a blue (dimers, top) or gray 
(spins, bottom) scale, respectively. 
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(a) ferromagnetic 
(t = -0.5) 
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FIG. 7. Direct comparison of dimer correlations - in an en- 
larged illustration of the top left hexagon of Figs. 6(b) and (c) 
- shows a small relative enhancement of correlations on two 
bonds for the disordered phase (V — 0.45 and original lattice 
L — 12). The upper bond is again the reference bond for all 
dimer correlations. 



range order. 36 ' 37 To check for this kind of order we calcu- 
lated higher-order correlation functions of the spin vari- 
ables, i.e., fourth-order correlations using improved esti- 
mators in the QMC simulations: 38 



(SjSjSfcS;) 



(SiSj-XSfcS,), 



(17) 



where i and j index sites at one nearest-neighbor bond 
and k and I at another nearest-neighbor bond. 

Figure 6 shows these correlations on a representative 
lattice where the strength of the correlation is given in 
a color code (blue scale, top) and the distance of the 
two bonds i-j and k-l is given by the distance of each 
bond to the top left reference bond. In addition, the S z 
correlation functions are shown in a gray scale (bottom 
scale) on the sites also with respect to the top left site 
of the lattice. In Fig. 6(a) the values of these correla- 
tions are shown for a Neel-ordered state where the S z 
correlations oscillate for different sublattices and show a 
constant nearly maximal amplitude. The dimer correla- 
tions are small and show no sign of ordering. For param- 
eters inside the in-plane ferromagnetic region [Fig. 6(b)] 
neither in the S z spin nor in the dimer correlations can 
any signature be detected, i.e., spin correlations drop 
rapidly to zero and dimer correlations adopt a constant 
distance-independent value. The same applies for the 
disordered region [Fig. 6(c)] and only a minor detail dis- 
tinguishes the two calculations: Apart from the differ- 
ent scales (cf. numbers at the upper scale) , which are ex- 
plained by the different strength of the quantum fluctu- 
ations (t — —0.1, — 0.5), we identify a small difference by 
comparing the relative values of the dimer correlations 
inside the top left hexagon shown enlarged in Fig. 7. In 
particular a slight enhancement of the dimer correlations 
on the bonds neighboring the opposite bond of the ref- 
erence bond (top) compared to the correlation on the 
opposite bond itself is seen. Thus an extremely short- 
ranged ordering of dimers is observed in the disordered 
phase, which is absent in the ferromagnetic state. 

Since we have found no finite signal for any magnetic 
or dimer order parameter we interpret the finite region 
in the phase diagram as disordered. To learn more about 
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FIG. 8. Spectrum of energy gaps AZ?k = — Eq for the 
Stotai = subspace for N = 34 at V = 0.45 and t = -0.1. 
The smallest gap at k = is stable in a finite-size analysis. 



the character of this phase, an investigation of the low- 
energy spectrum is necessary. Therefore, we performed 
an ED of the lattice model for a set of parameters in the 
disordered region (V = 0.45 and t = — 0.1). In Fig. 8 we 
show the energy differences A£ k = E k — Eq in the S* otai = 
subspace. This measures the gap from the ground-state 
energy Eq (which belongs to the k = subspace and its 
scaling is shown in the top panel of Fig. 9) to the lowest 
eigenvalues E^ in the different k / subspaces or to the 
first excited state in the k = subspace. The smallest gap 
in this same subspace stays finite according to a finite- 
size analysis shown in Fig. 9 (middle, AF gap). 



^total 



We also calculated the lowest eigenvalues for higher 
subspaces and found that the gap between E and 
the lowest eigenvalue in the S* otal = 1 subspace vanishes 
in a finite-size scaling analysis for system sizes N = 18 - 
34 (plotted in the bottom panel of Fig. 9). This indicates 
a ferromagnetic correlation, which can also be seen in the 
correlation functions of the S xy components in the ED. 
However, from the scaling of the ferromagnetic gap it is 
obvious that the ED suffers from finite-size problems for 
the small lattices, which were investigated since we could 
exclude ferromagnetic order from the QMC simulations. 



In spite of finite-size effects, the clearly finite gap in the 
S z otal = subspace found by ED in the disordered phase 
for the present model allows us to exclude a topologically 
ordered state. The reason is that, on a two-dimensional 
periodic lattice with an aspect ratio close to one (i.e., on 
a torus with similar circumferences in both directions), 
a topologically ordered state would be accompanied by a 
fourfold ground-state degeneracy. 39 
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FIG. 9. Scaling of the ground-state energy and two gaps with 
the system size N from ED at V = 0.45 and t = -0.1. The 
antiferromagnetic gap (AF, middle) inside the St ota x = sub- 
space stays finite; however, the ferromagnetic (FE, bottom) 
gap between the ground-state energy and the lowest eigen- 
value of the 5 t z otal = 1 subspace scales to zero. 



ters we calculated higher-order correlation functions and 
excluded any long-range order of dimer configurations. 
This absence of any finite order parameter is in agreement 
with earlier calculations for the isotropic model, 8,10,11 al- 
though the stability region of this disordered phase differs 
at a quantitative level. A further investigation of the low- 
energy spectrum using ED was hampered by finite-size 
effects but yielded a finite antiferromagnetic gap, which 
indicates absence the of topological order. For isotropic 
models the disordered state with a finite gap would be 
referred to as a gapped spin-liquid phase. In principle, 
similar finite-size effects as in ED are also present in the 
QMC simulations, but here they can be overcome by cal- 
culating larger systems. This underlines again the utility 
of introducing ferromagnetic quantum fluctuations. 

The phase diagram is very similar to the one 
found for the square lattice with nearest- and next- 
nearest-neighbor interactions, which was investigated for 
anisotropic exchange terms in Ref. 15. However, the sta- 
bility region of the antiferromagnetic phases is reduced on 
the honeycomb lattice, which is explained by the lower 
coordination number of the lattice, thus the influence 
of quantum fluctuations is enhanced. In particular, the 
critical value of the ratios t r /V r obtained by QMC simu- 
lations at which the antiferromagnetic phase boundaries 
split and the disordered phase emerges is smaller in the 
present case. Furthermore, the direct transition line be- 
tween the two antiferromagnetic states calculated via SE 
shows a steeper slope for the square lattice which also 
hints at a larger stability of this direct transition in that 



case 



15 



V. CONCLUSION 



We presented the phase diagram for an anisotropic 
Heisenberg model on the honeycomb lattice with up to 
third-nearest-neighbor interactions. The quantum fluc- 
tuations in the present model are chosen to be ferromag- 
netic in order to avoid the sign problem in QMC simula- 
tions. Since apart from the interaction J2 this sign can 
be absorbed by a sublattice rotation, one may hope to 
nevertheless gain insight into the isotropic model. 

The frustrating next-nearest-neighbor interactions 
suppress the Neel state, which is favored by the anti- 
ferromagnetic S z interactions between nearest and third- 
nearest neighbors. Thus a collinear state arises for large 
interactions between the sites on the same sublattice. 
The schematic transition lines were calculated by means 
of LSWs and the direct transition between the competing 
antiferromagnetic states was found to survive for small 
fluctuations by means of SE. In addition, we also per- 
formed QMC simulations to determine the phase dia- 
gram and in particular to gain insight into a region of 
the phase diagram without any conventional order. Since 
an earlier ED study of the isotropic model reported the 
appearance of valence bond crystals 10 for certain parame- 
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